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Abstract 

We address the problem of navigating a set (fleet) of aircraft in an aerial 
route network so as to bring each aircraft to its destination at a specified 
time and with minimal distance separation assured between all aircraft at all 
times. The speed range, initial position, required destination, and required 
time of arrival at destination for each aircraft are assumed provided. 

Each aircraft’s movement is governed by a controlled differential equation 
(state equation). The problem consists in choosing for each aircraft a path 
in the route network and a control strategy so as to meet the constraints 
and reach the destination at the required time. The main contribution of the 
paper is a model that allows to recast this problem as a decoupled collec- 
tion of problems in classical optimal control and is easily generalized to the 
case when inertia cannot be neglected. Some qualitative insight into solu- 
tion behavior is obtained using the Pontryagin Maximum Principle. Sample 
numerical solutions are computed using a numerical optimal control solver. 

The proposed model is first step toward increasing the fidelity of continuous- 
time control models of air traffic in a terminal airspace. The Pontryagin 
Maximum Principle implies the polygonal shape of those portions of the 
state trajectories away from those states in which one or more aircraft pair 
are at minimal separation. The model also confirms the intuition that, the 
narrower the allowed speed ranges of the aircraft, the smaller the space of 
optimal solutions, and that an instance of the optimal control problem may 
not have a solution at all (i.e., no control strategy that meets the separation 
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requirement and other constraints). 
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1. Introduction 


Problems in coordinated motion planning for multiple agents can be 
roughly classified into two disjoint categories, decoupled coordination (each 
agent’s motion is planned separately, then the plans are reconciled), and 
centralized coordination (all the agents’ motions are planned simultaneously, 


with the interaction constraints considered from the start) ( Ghrist et ah 
2005). The problems considered in this paper fall in the latter category. 


Centralized coordination of multiple agent motion has been approached us- 
ing various types of mathematical models, discrete (see, for example. Luna 


2011 


2006 


and references therein) and continuous (see , for example, Hu et al. 20031. 


Jung and Ghristl 120081: iGhrist et al.l 120051. and references therein) . A 


review of research on multi-robot coordination problems can be found in 
Ghrist et al. ( 2005 . section 1.1). 

In a number of coordination problems, the moving agents are confined 
to a transportation network (also known as roadmap coordination space ; see 
Ghrist et al.l 2005 ). A general mathematical model of a transportation net- 
work is a multigraph ( Papadimitriou and Steiglitzl 1998 1 with the vertices 
being points in a Euclidean space and edges being parametrized curves con- 
necting pairs of vertices. Examples of such networks include railroad net- 
works for trains, railroad networks for industrial robots, trolley and tram car 
networks, and airspaces with fixed nominal routes. 

A subclass of network-confined coordination problems consists of those 
where the agents’ paths are not given but sought as part of solving the prob- 
lem. In such problems, the system exhibits behaviors both continuous (the 
agent’s motion along an edge) and discrete (an agent’s choice between two 
edges emanating from the same vertex). This coupling of the two behaviors 
suggests hybrid control systems (HCS) (jSussmannl . 1999) as a suitable class 
of models for approaching the problem. Hybrid systems have been applie d to 
various problems o f transportation (e.g., highway traffic: iMunoz et al.ll2 003: 


Baven et al.l I2004T). and i n particular to aerospace problems (ISastrv et al 


1995 : Hwang et al. . 2006 : Baven and Tomlin . 2003h . 
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Although the HCS model proposed herein is in principle applicable to 
all types of moving agents whose motion is to be coordinated in a given 
transportation network, this paper will be held to the specific context of Air 
Traffic Management (ATM). Thus, the moving agents here will be aircraft, 
and the central problem will consist of finding routes and speed advisories for 
a set of aircraft, subject to initial conditions, separation constraints, speed 
range constraints, and an arrival schedule. The role of the transportation 
network will be played in this context by an aerial route network (for brevity, 
it will be called simply route network ), i.e. a collection of pre-defined aerial 
routes, to which all aircraft motion is to be confined, in a given airspace. 

This paragraph and the next two describe the context for the HCS prob- 
lem central to this paper. The context here is a terminal airspace that consti- 
tutes a relatively small part of the U.S. National Airspace System (U.S. NAS). 
A typical setting is a Terminal Radar Approach Control Facility (TRACON) 
operation that spans 0-2 hours. In this context, the role and, therefore, com- 
putational demands, of the routing problem (i.e., the combinatorial problem 
of finding for each aircraft a walk in the multigraph that represents the 
route network) depend on the current ATO practices, on the prognoses for 
NextGen, and on the density (see below) of the multigraph. These depen- 
dencies are s ummarized in Ta_ble [T] Henceforth, the term dense refers, as in 
dense graph ( Corrnen et all 1199011 . to a high ratio between the number of 
edges and the squared number of vertices; sparse is the opposite of dense. 

In today’s operations, the situation is as follows. Since these operations 
rely on air traffic controllers and pilots to keep the collision risk between air- 
craft below an acceptable threshold level, a sparse route network in terminal 
airspace aids the situational awareness required for controllers and pilots to 
perform separation tasks. The limited number of arrival routes in a given 
airport runway configuration allows to route aircraft hours in advance of ar- 
rival into the terminal airspace, thus removing the real-time urgency from 
the routing problem. While air traffic controllers instruct aircraft to turn 
away from the assigned route for separation from other aircraft, the general 
flow of aircraft (and situational awareness) is retained. 

Whether future ATO as proposed by NextGen and SESAR use a sparse 
terminal route network, as they do today, or a denser route network to achieve 
higher runway utilization, is to be determined. In general, a denser route 
network may increase runway utilization and lessen the inefficiencies that 
would have arisen from handling operational uncertainty with limited routing 
options. Achieving this at the expense of safety, however, is an unacceptable 
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Table 1: Implications of the complexity of the routing problem. 


trade. 

For the HCS problem central to this paper, one is given the following 
data: 

1. A directed multigraph G = ( V, E ), with vertex set V and edge set E, 
each vertex v G V being a point in a Euclidean space E of dimension 
2 or 3. If e € E is an edge from vertex V\ to vertex u 2 , then the nom- 
inal route segment from waypoint V\ to waypoint v 2 is a curve in E, 
connecting u-i t o v 2 . All such curves will henceforth be assumed rec- 
tifiable ([Shilov, 1996) and capable of a parameterization which is con- 
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tinuous and piecewise continuously differentiable. A cusp in the curve 
can be traversed with the assumption (made throughout this paper, 
but capable of relaxation) that inertia is neglected, and approximately 
smoothed if inertia i s to be taken into account. A graph-theoretic path 
f Papadimitriou and Steiglitz) . 1998 1 in G is, therefore, associated (and, 
henceforth, identified) with a spatial path th at ca n be traversed by an 
aircraft. A vertex of G of indegree > 2 flPapadimitrion and Steiglitz . 


1998 ) (resp., outdegree > 2) corresponds to two or more route segments 


4. 


merging (resp., diverging). The modeling framework below imposes no 
restrictions on the outdegree or indegree of a vertex. 

A finite set 

A — A} 

of aircraft a G A in G. If aircraft a is moving along a path in G, then 
the aircraft’s position is specified by the arc length coordinate x a along 
the path. 

For each aircraft a G A, a specification of its initial position x INIT ' ,a , 
required destination x DEST ’ a , and the required time t DEST ' a of arriving 
at the destination. Here x INIT ' ,a and x DEST,a are points in G, each 
point specified, for example, by an edge in G and a fractional distance 
along that edge. 

The inertia- free state equations ( Korn and Korni 196lh (henceforth the 
dot denotes differentiation with respect to physical time t) 


x = s 


CL G A 1 


where the s“’s are the corresponding speeds, describing the motion of 
those aircraft a that have not yet reached their destination. In what 
follows, and with the details provided below, the coordinates x a will 
play the role of state variables ; the speeds s“, of the control variables. 

5. State constraints: the separation requirement for each pair of aircraft. 
This requirement is described mathematically, in terms of the coordi- 
nates x a , in section [3j 

6. Control constraints: bounds on the speeds s“. 

7. A cost functional , specified below. 

The problem, defined in detail below (definition 14. lj) as the Scheduled Routing 
Problem , consists in finding for each aircraft a G A a path p(a) in G from 
x iNiT;a x DEST;a anc [ a con f ro i strategy Sp( Q )(f), i.e. a function from the time 
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domain of the model to the set of admissible controls, such that the resulting 
movements x a (t) along the corresponding paths constitute a state trajectory 
that satisfies the above state and control constraints and that minimizes the 
cost. 

In this paper, we use an HCS framework to formulate a model specialized 
to the above problem. The main contributions of this model are as follows: 


Reduction of the problem to a special case of an HCS where each so- 
lution trajectory lies in only one control mode. 


A clear application of Depth-First Search ( Cormen et all 1990) to 
search through the control modes as economically as possible, given 
the possibly exponential size of the problem. In the worst case when 
every aircraft can be assigned to any of the paths, and when each edge 
can serve as a path by itself, the number of aircraft-to-path assign- 
ments, i.e. of functions /i : A — > E, is |£)I‘ A L The latter quantity, 
however, is a crude overestimate: the set of an aircraft’s paths from 
initial position to destination is unlikely to involve a large subset of E. 
Nevertheless, if every aircraft has even as few as two such paths, then 
the number of routings is 2^1, which is still exponential in the number 
of aircraft. This exponential behavior is a challenge and is mitigated 
by two factors: 


— the applicability of parallel computing, and 

— the operational constraints of ATM (at the time of this writing) 
that restrict the number of possible aircraft-to-path assignments. 

These factors are discussed briefly in section 14.51 


Reduction of each control mode to a problem in classical deterministic 
optimal control, which allows, at least in principle, application of the 


funda mental results of Pontryagin (IBoltvanskii et all Il962l) and Bell- 


man f Bellmanl . 2003 1. and of the numerica l algorithms t hat have been 
developed and implemented (Pytlakj, 19991 : Fab ienl. 2011). 


• A natural way to capture an aircraft’s exiting the system; see Remark 
HU below. 


These contributions together allow for parallel computation of solutions: the 
classical optimal control problems corresponding to different control modes 
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0 0 


can be solved in parallel, and their obtained minimal costs values compared. 
Furthermore, the Depth-First Search algorithm itself admits a parallel im- 
plementation ( Rao and Kumar . 19871) . 

The hybrid model is formulated in section [2] The classical deterministic 
optimal control corresponding to a given control mode is formulated in section 
Numerical solutions to some instances of the problem are given in section 


2. An HCS formulation 

The HCS defined in this section will be instrumental in a precise formu- 
lation of the Scheduled Routing Problem. Assume the data -|7J), listed in 
the 3rd paragraph of section [Q 

• For each aircraft a G A, let e /7V7T;a denote the edge occupied initially 
by aircraft a. 

Remark 2.1. However, since a given routing assigns each aircraft to 
a specific route, which is a parameterization of a path in G ( see sec- 
tions \4 ■ S\ for more detail) and to which the aircraft’s movement is 
restricted, the only way for one aircraft to overtake another is if their 
two routes permit. In particular, an aircraft may not overtake another 
on the same edge. 

• Let V (e INIT ' a , e DEST ’ a ^ be the set of all paths in the multigraph that 
begin with the edge e INIT ’ a and end with the edge e DEST ’ a that contains 

nDEST ;a 
X 

Definition 2.1. The length of a path p G V (e INIT ' a , e DEST > a 'j W M be 
denoted l(p) . 

• Define a control mode p as a mapping that assigns each aircraft a to a 
path in V (e INIT,a , e DEST '' a \ . In more detail, fi is a mapping from the 
set A of moving aircraft to the union U a V (e INIT ' a , e DEST;a ) such that 

p,(a) G V [e INIT ' a , e DEST ’ a ) for each a G A 
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• For aircraft a , each path /z(o) G V (e INIT '°‘, e DEST ’ a ) is parameterized 
by arc length. For computational convenience, the arc length coordi- 
nate increases along the path, with the destination coordinate x DEST ' ,a 
being zero for each a. Thus, x" G [—/(//(a)), 0], and 

x DEST ’ a = 0 

This convention ensures that an aircraft’s destination is a vertex in the 
route network and, furthermore, that it is the same vertex in all control 
modes. 

• Each aircraft a in each control mode p. is required to reach its desti- 
nation x“ = x DEST ' a = 0 at a prescribed time t DEST ' ,a . Upon reaching 
destination, the aircraft is no longer in the model; this is reflected in 
the restriction on the time domains of the individual state equations 
©, stated below. 

• In each /z, have the arc length coordinate x“ evolve according to the 
state equations 

x“(t) = slit) for 0 < t < t DEST ’ a , aeA , (1) 

where is the control variable corresponding to the aircraft’s speed 
of motion along the path. In this formulation, the state equations are 
imposed over different (albeit overlapping) time domains. The problem, 
however, will be converted below to one where all state equations are 
imposed over the same (rescaled) time domain. 

• For each /z and each a , impose the arrival requirement 

DEST;ct 

• For each /z and each a , impose the speed ranges 

sf IN ’ a < s“ < sf? AX '' a (3) 

3. The geometry of separation constraints 

In some transportation types, including aircraft and trains, every pair 
of moving agents must be-as a safety measure-separated by a distance 



x 


DEST\a 


( 2 ) 


no smaller than a pre-determined minimal separation. For aircraft in 
the U.S. Air space, the minimal separation requirements are defined by 
Order JO 7110. 65U (2012) and depend on numerous factors, including 


airspace type, air traffic automation system s in use, and aircraft weight 
classed (the classes defined by Order JO 7110.65U ( 2012 1 are: Small, Large, 
Heavy, B757). The separation requirement will be a key constraint on the 
state variables in the Scheduled Routing Problem, formulated below. Recall 
(remark 12. lj) that the movement of an aircraft is restricted to the assigned 
route. 

No attempt is made in this paper to capture all such requirements in 
detail (see, however, section 16.31 for a discussion of location dependence in 
separation requirements). Instead, we will use conservative approximations, 
addressing only the following asymmetry: if two moving aircraft are in-trail 
(i.e., one is directly following the other along a route segment which is not 
necessarily in a horizontal plane), then the minimal separation can depend on 
the weight class of the leading and trailing aircraft. To capture this potential 
asymmetry, for each pair aq,a 2 of aircraft with the first one leading, we 
introduce the minimal separation r Q1 . Q2 . If the asymmetry takes place, it can 
be written 


Q!i,Q2 


r, 


(* 2 , 0:1 


( 4 ) 


We now calculate the set of all the “forbidden” states, in a control mode 
p of a hybrid system described above; i.e., the set of all states such that 
at least two aircraft would have violated the separation requirement if the 
system were to enter that state. The scenario shown in Figure [1]A has two 
aircraft on two different rectilinear edges, which need not lie in a horizontal 
plane, with a common vertex and no specified orientation. (If the edges are 
curvilinear with low curvature near a common vertex or intersection, these 
portions can be approximated by linear segments; otherwise, the analysis 
becomes considerably more complicated.) 


Remark 3.1. Since edge orientation is not specified, Figured describes four 
cases: both aircraft are moving toward the common verted, both moving away 


2 The term weight class is used here because it is an FAA term lOrder JO 7 1 10.65U 


( 201 2 1 . 

3 In ATM, such a vertex is called a merge point. A merge point either must have a t hird 
edge emanating from it (i.e., must have outdegree ( Panadimitriou and Steiglitzl . 1998) at 
least 1), or must be a final point of arrival for both aircraft. 
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(A) 


(B) 


Figure 1: Aircraft 1,2 on their respective rectilinear edges ei,e 2 , which share a common 
vertex, taken as the origin 0 in R 2 . The orientation of the edges is not specified. (A) 
The unit vectors ai,a 2 are collinear with the respective edges, but their directions do not 
necessary agree with the edges’ orientations. (B) With suitably chosen scalars coefficients 
Ci, C 2 , the vectors and C 2 a 2 are the respective position vectors of the two aircraft. 


from the common vertex, and two more cases in which one aircraft moves 
toward, and the other away from, the common vertex. 

The case when the two aircraft are on edges that do not share a common 
vertex yet are sufficiently close together ( possibly intersecting) to allow vio- 
lation of the separation requirement is treated in the last paragraph of this 
section. 

We will use the Euclidean inner product (•, •) and the corresponding norm 
1 1 • 1 1 in the 2-D space containing the two edges. Pick the coordinate system^ 
with the common vertex as the origin and with the unit vectors ai,a 2 as 
the basis vectors that, regardless of the edge orientations , point from the 
origin toward the respective aircraft. With suitable scalars ci, c 2 , the vectors 
Cia! and c 2 a 2 are the respective position edges of the aircraft. The squared 
distance between the two aircraft is 

||ciai - c 2 a 2 1 1 2 = (ci) 2 + (c 2 ) 2 - 2cic 2 (ai,a 2 ) (5) 

Equating the latter expression to the squared minimal separation, say, r 2 2 , 


4 This coordinate system is chosen here because it results in the simplest form of the 
squared distance between the two aircraft (the form we obtain here is the left-hand side 
of©). Other systems can be used, possibly resulting in a more cumbersome algebra. 
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we obtain the equation 

( c i) 2 + ( c 2 ) 2 — 2cic 2 (ai, a 2 ) = r 2 2 (6) 

of an ellipse in the CiC 2 -plane. The corresponding set of conflicting sets is 
described by the elliptical sector obtained by intersecting the ellipse-bound 
region 

(ci) 2 + (c 2 ) 2 - 2cic 2 (ai, a 2 ) < r 2 2 

with the open octant c 2 > C\ > 0, corresponding to the case when aircraft 1 
is the one closer to the origin. In the other case (aircraft 2 is closer to the 
origin), the corresponding elliptical sector is obtained by intersecting 

( c i) 2 + ( c 2 ) 2 — 2cic 2 (ai, a 2 ) < r 21 

with the octant Ci > c 2 > 0. The role of the angle 9 between the edges 
ei,e 2 in both sectors is the equality (ai,a 2 ) = 008(0). An example of two 
such sectors is shown in Figure [2l The asymmetry of the gray-shaded region 


Figure 2: An example of two elliptical sectors in the ciC 2 -plane corresponding to conflicting 
states. 

about the dashed diagonal is the asymmetry (JH). 

In each of the four cases listed in Remark 13.11 the respective continu- 
ous state coordinates x^x 2 ^ of aircraft 1,2 in control mode /i map to the 
coefficients Ci, c 2 , as follows (recall the notation introduced in definition 12. lft : 

1. If both aircraft are moving toward the common vertex, then = 
/(e Q ) — c Q for a — 1,2. 

2. If both aircraft are moving away the common vertex, then = c a for 

a = 1 , 2 . 
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3. If aircraft 1 is approaching, and aircraft 2 going away from, the common 
vertex, then x 1 = l(e i) — ci, x 2 — c 2 . 

4. If aircraft 2 is approaching, and aircraft 1 going away from, the common 
vertex, then x 1 — ci, x 2 ^ — l(e 2 ) — c 2 . 

If 9 > 90°, then in the last two cases /i allows only one in-trail sequence, so 
the minimal separations used for the two sectors in Figure [2] are equal. If the 
two aircraft 1, 2 are on one and the same edge, then the set in the CiC 2 -plane 
of the conflicting states appears as in Figure [3] (the asymmetry about the 
dashed diagonal corresponds to (JH). The mapping from the continuous state 



Figure 3: An example of two stripes in the ciC 2 -plane corresponding to conflicting states 
of two aircraft on the same edge. 


coordinates x^x 2 to the coefficients ci,c 2 is constructed analogously to the 
above four cases. 

The above calculation is illustrated, for an example of two aircraft, in 
FigureJdJ Each di screte mode’s set of conflicting states is shown as a connected 
(I A dams and FranzosaU2008l) gray region. For dimension A above 2, one must 
compute for each pair of aircraft the set of states violating the separation 
requirements. Each such set is a cylinder, or union of cylinders, with the 
base shaped as shown in Figure HP, in the total state space We 

note that the set o f all separa tion-violating states in U M X (J is cylindrical in 
the sense of ( .Tung and Ghristl . 2008, Definition 2.2), the latter definition a 
key requirement for the applicability of a number of theoretical results of 


Jung and Ghrist (20081). 


A case that remains to be treated is that of two aircraft, labeled here 
1 and 2, on edges that do not share a vertex, yet allow the two aircraft to 
violate the separation requirement. The two edges may or may not intersect. 
On denoting the position vectors of the edges’ initial vertices as hi,h 2 , the 
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Figure 4: An example of two aircraft whose paths in the route network are prescribed and 
overlap. The black star shows the beginning of the overlap in (A) and the corresponding 
state (both aircraft being at that point) in (B); the white star, the end of the overlap in (A) 
and the corresponding state (both aircraft being at that point) in (B). The system, shown 
in (A), has 7 discrete modes (B) with both aircraft in the rou te networ k. Each mod e’s set 
of separation- violating states, shown in (B) as a connected (I Adams and Franzosa . 2008 ) 
gray region, is “glued” to some of the others. The result of the gluing is the connected 
region shown in (C). 


above calculation can be generalized to this latter case by calculating the 
squared distance between the aircraft using, instead of the left-hand side of 
(J5]) , the expression 

1 1 (hi + Ciai) — (h 2 + C2a 2 ) 1 1 (7) 

and updating the subsequent calculations accordingly. Note that taking both 
h[ , h 2 to be at the origin would reduce the general formula (J7J) to the special 
case depicted in Figure [D 

4. The Scheduled Routing Problem and the equivalent Stacked 
Scheduled Routing Problem 

For ease of exposition, we precede the general formulation of the problem, 
suitable for an arbitrary number of moving aircraft, with a specific, two- 
aircraft, example. 

4-1. An instance of the scheduled routing problem for two aircraft 

The initial state of a two-aircraft scheduled routing is shown in Figure 
[5jA; the required destinations (with possibly different required arrival times) 
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for the two aircraft, in Figure 03. In this Figure, 


a — 1 : e INIT ' 1 = ei, e DEST ;1 = e 8 , 
a — 2 : e INIT ;2 = e 2 , e DEST;2 = e 7 . 

The only paths that take the aircraft from initial location to destination 




(A) 


(B) 


Figure 5: An initial state (A) and required destinations (B) of a 2-aircraft set. 


(aircraft 1 from e\ to e 8 , and aircraft 2 from e 2 to e 7 ) are 


Pi '■ ei, e3, e6, e 8 ; P2 ; ei,e 4,e 8 ; P3 : e 2 , e 3, e 6, e 7 ; Ah : e 2 , e4, e 7 . (8) 


Aircraft 1 can take either pi or p 2 ; aircraft 2, either p 3 or p 4 . Consequently, 


V{ 


’ INm \ e DEST ’ 1 ) = {p u p 2 }, V {e INIT ’\ e DEST ’ 2 ) = {p 3 ,p,}. 


No other paths can be taken. The paths Pi, ■ ■ ■ ,P 4 are shown, in that order, 
in Figure EO We obtain the control modes ... ,^4, defined as follows: 


a 

Ah (a) 

Ah(a) 

p 3 (a) 

Ah (a) 

1 

Pi 

Pi 

P2 

P2 

2 

P3 

Pa 

P3 

Pa 


The state space corresponding to each control mode p is a rectangle consisting 
of those state vectors (x^, xfy which are compliant with the arc length bounds 
and separation constraints. 

The above system is subject to the operational requirement (J2]), here for 
a G A = {1,2}, that aircraft 1 and 2 arrive at their destinations at times 
j-DRST: 1 , j-DEST]2 ^ reS p ec tively. 
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Figure 6: The four paths ([5)1 . in that order, for the problem in Figured] 


In each of /q , u 4, the pairwise conflict zone is simply connected 
( Adams and Fra/nzosa . 2008h : i.e., consists of only one connected component 
and has no “holes” . Consequently, once each of the conflict zones, of the form 
shown in Figure [4p, is approximated by an ellipse-bounded region, the state 
spaces for /ii, /i 4 have the topology shown in Figure [7]A.. In each of p 2 , /i 3 , the 
paths of the two aircraft have two crossings (regarded here as short overlaps), 
hence the pairwise conflict zone has two connected components (Figure [7)3). 


Each control mode p is subject to the initial condition 

x"(0) = x INIT,a for a G A (9) 

The simplifying assumption underlying ([9]) is that both aircraft start their 
movement simultaneously. This assumption can be relaxed and is used here 
for mathematical simplicity only. 

Finally, for each control mode p, we can specify a cost functional that suits 
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A 
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6 - 


1 




(A) (B) 

Figure 7: Topology of approximated conflict zone in the state spaces (A) for /zi and /i 4 
and (B) for and H3- The black dot (near top right), the required destination coordinate 
pair (x EEST ’\x EEST ’ 2 ). 


the goals implied by the context of the specific application. One example, 
impractical but chosen for simplicity, is a cost functional that equals the 
average (squared) speed of the aircraft during the flight: 

t DEST-,a 

E / Of* < 10 > 

aeA Jo 

Thus, for each control mode p, we have an optimal control problem with 
state equations subject to the initial condition ([9]), the control constraints 
that specify permissible value ranges for the speeds s“, and the following 
additional constraints: 


The arrival requirement (j2j). 

The separation requirement, defined by constructing a function 

9° EPm ’ a2 (KUA) 

of the stat^f] such that the pairwise conflict zone (shown as a gray- 
shaded region in the appropriate panel of Figure [7]) is exactly the set 


5 This construction has already been carried out: the latter function is to be defined as 
the squared distance 0, where all the quantities with subscript 1 correspond to cu; all 
those with subscript 2, to < 22 . 
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of states satisfying the functions 


SEP,ai,a2 \ r 2 

— ' ai,Q:2 * 


(ii) 


For the pair (an = 1, a 2 = 2) (and, in the general scheduled routing problem, 
for every pair (ai,a 2 )) of aircraft, inequality (ITT]) defines a union (denoted 
X“ 1;Q2 ) 0 f regions in the state space of p, each region bounded by an ellipse- 
cylindrical (“tube-shaped”) hypersurface. The separation requirement thus 
translates into the requirement that a solution (the state trajectory) is dis- 
joint from the interior of every such region for every pair of aircraft. 

4-2. The Scheduled Routing Problem: a general formulation 

The central problem of this paper, which will be called the Scheduled 
Routing Problem , can now be stated as follows. 

Definition 4.1. (Scheduled Routing Problem) Given a set A of aircraft 
moving on a route network G = (V, E) subject, in each control mode p, to 
the state equations IfTf), the initial condition (GJ), the state constraints / fill) . 
and the control constraints fS j), find 

(a) a control mode p and 

(b) a corresponding control strategy (Sq(t))q, such that the corresponding 
state trajectory x M (f) = (xf(t)) a£ ^ 4 satisfies (dj) and minimizes the to- 
tal cost over all the control modes. 

Part (jaj) of the Scheduled Routing Problem consists of routing the set A 
of aircraft; i.e. , of choosing for each aircraft a path from origin to destination. 
A solution found for part ((a]) heavily affects the process of solving part (lb|) : 
there may be control modes p for which part (JbJ has no feasible solution. 
Some ways to approach part (Jaj) and the corresponding computational costs 
are discussed in section 14.51 

Part (0 of the Scheduled Routing Problem is an optimal control problem 
with a Lagrange cost function. The problem, however, has two non-standard 
features that hamper application of classical optimal control theory and nu- 
merical computation of solutions. One feature is the presence of intermediate 
constraints: in the arrival requirement ([2]), all time instants t DEST ' ,a except 
the latest one are interior points in the time domain of the problem. The 
other is the following non- autonomous behavior: an aircraft, once at desti- 
nation, no longer “participates” in the constraints or in the cost. We now 
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use a formalism similar to that in IDmitnik and Kag anovich (2 01111 to reduce 
this problem to a classical, optimal control problem, which is autonomous if 
its cost functional is. In the rest of this section, the subscript p is dropped 
for brevity. 


4-3. The Stacked Scheduled Routing Problem ( SSRP ), equivalent to the 
Scheduled Routing Problem 

Let (a q )^ =1 be an ordering of the aircraft by their arrival times (arranged 
in nondecreasing order). Let t 0 = 0, and for each q > 0 let 

j. / .DEST\a q 


Define a new (normalized) time r G [0,1] and, for each of the intervals 
[' t q , t q . |_i], q — 0, . . . , A — 1, introduce the following: 

• the state variable p 9 (r), which will play the role of “physical time” in 
the interval [t q ,t q+ i\: 

tq 5; Pqi' 1 ') ^ t q + 1 , 

• the state variables y q {r ), aed, related to the above x a (t ) by 

Vq(r) = x a (t ) if p q (r) = f; (12) 

• the vector notation y 9 (r) = (y q {r)) a&A \ 

• the control variables s q (r), related to the above s a (t ) by 

Sq(r) = s a (t ) if p q (r) = t , (13) 

• the control variables z q > 0, which represent the “rate of flow of physical 
time” with respect to the normalized time r. 


The above definition of the Lagrange cost reflects the assumption that an 
aircraft, once at destination, “disappears” from the system, in the sense of 
being no longer subject to the separation requirement with the other aircraft. 
From ([T]), (TT2]h and (fT3jh one readily obtains the state equations 


£yq q ' = z q Sq ql for q' > q + 1 
!SPq = z i 


OL G wT. 


(14) 


18 


Finally, the initial conditions and arrival requirement, together with the re- 
quirement that physical time and state change continuously when passing 
from one interval [t q ,t q + 1 ] to the next, translate to the endpoint constraints 


ho(0) 

X INIT;a ( see ^Jgj^ 

(a) ' 



Pqi 1 ) 

Pg+l(0) 

(b) 



Pq+l(0) = 

t q +l 

(c) 

> 0 < q < A 

(15) 

y?'( 1) = 

l/qii (°) for q' > q + 1 

(d) 



i) = 

^,DEST;a q +i 

(e) , 



Conditions (flSlbc) ensure the continuity of physical time flow; 

conditions 

flUld), the continuity of 

an aircraft’s motion. 





Remark 4.1. The state equations [Lfy and the endpoint constraints [TRd) 
come with the restriction q' > q + 1 because they are imposed, in agreement 
with m, only for those aircraft that have not yet reached their destination: 
by the end of the q-th time period, the first q aircraft have reached destina- 
tion. Thus, upon reaching destination, an aircraft is excluded from the model, 
and is no longer represented by a state equation or subject to separation con- 
straints with the other aircraft. This eliminates the necessity to continue 
modeling an aircraft whose role in the model has already been fulfilled. 


Throughout the rest of this paper, the choice (TTOT) of the cost functional is 
assumed in all the numerical examples. Other choices of the cost functional 
are discussed in section |6j Since the time intervals between each consecutive 
pair of arrivals are modeled, in this latter formulation, as if they were occur- 
ring simultaneously (“stacked” upon one another), the following definition 
will be adopted. 


Definition 4.2. (a) To each control mode p, associate the newly obtained 
optimal control problem consisting of: the state equations 0. endpoint 
constraints as , cost functional as. the separation constraints on the 
variables y q q , the control constraints corresponding to (0) specifying the 

Oi / 

speed ranges on the variables s q q , and the positivity constraints 


z q > 0, 

will be called a /r-stacked optimal control problem. 
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(b) The set of all fa-stacked optimal control problems will be called a Stacked 
Scheduled Routing Problem (SSRP). 

(c) Of all the optimal solutions to all the fi-stacked optimal control problems, 
a solution achievinq a lowest cost is called an optimal solution to the 
SSRP. 


Thus, an SSRP consists of a collection of optimal control problems, and an 
optimal solution to the SSRP tells not only how quickly the aircraft are to 
move, but also how they should be routed. The SSRP is equivalent to the 
Scheduled Routing Problem (definition 14. ip . 


4-4- Implications of the Pontryagin Maximum Principle for the SSRP 

The assumption that we are in a specific control mode fi is still in force. 
Denote by f Vq . f Pq the costate variable cl for y q , p q . In those states where none 
of the state constraints is active, the Hamiltonian for each control mode of 
the SSRP is 

H = ~fo + ^2 ZqSq q iy q + ^ Z q f Pq , 

q-,q’>q + 1 Q 

where / 0 is the performance index (running cost) corresponding to the cost 
functional (TT01) . Since H does not explicitly depend on any of the state vari- 
ables, it follows that the costate variables are constant along the trajectory 
portions clear of state constraints and, consequently, the maximization of H 
in each such state is a problem of static maximization. This, in turn, im- 
plies the existence of an optimal state trajectory in which these portions are 
segments of a straight line on each [t q ,t q +1 ], If, furthermore, there is only 
one optimal trajectory, then its portions away from obstacle boundary are 
necessarily rectilinear. 

With the state equations ([I]) and the set (J3J) of control constraints, the 
set of all states reachable from a given state y° = {y p a )a£A is the pointed 
polyhedral cone that consists of all states y p = (yf) a eA satisfying 


0 < (»“ 



a MIN-a i 

jV 

MAX;a 2 


< V 


Q2 



< 



c MAX-, ai 



MIN\d2 : 
Sy 


«1 7 ^ « 2 - 


6 Also referred to as adjoint v ariables in literature on Pontryagin’s Maximum Principle; 
e.s.. see {Korn and Korn. 1961. section 11. 8-2(a)). 
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This cone has vertex yjj and is the intersection of the half-spaces 


and 


a MAX ;ai 
^ /- 

y V ~ — ~MIN-,a 2 


a.2 ^.0;a2 


s MIN;a i 

,.“2 _ 7 ,0;O2 ■> _/f 

y/i y/i — MAX', a.2 


{y^ 1 - y°f ai ) > « 17 ^« 2 , 




The narrower the speed ranges (J3]) , the narrower the cone, and the smaller 
the portions of an optimal trajectory that lie on the boundary of the obstacle 


I | v a i’“ 2 

u ai?ia2 y '-(i 

This suggests that, with narrow speed ranges, optimal trajectories for the 
Stacked Scheduled Routing Problem can be well approximated by piecewise 
linear curves. 


4-5. An approach to computing solutions and an analysis of the associated 
costs 

The total number of the control modes /! possible at the outset (i.e., 
before checking for those /j not allowing a feasible control strategy) is 

| IN IT-, a e DEST-,a 

a&A 

Paths which are, or contain, cycles are allowed in the model and can be 
desirable in some applications, e.g. in air traffic models where aircraft may 
be sent into a holding pattern to absorb delay. With regard to the SSRP 
specifically, the following holds true: 

Remark 4.2. Every solution to the control problem described above will cor- 
respond to exactly one of the control modes fi. This removes two difficulties 
associated inherently with hybrid systems and absent from classical control 
systems, the risk of excessively frequent switchings of control mode, and the 
necessity for “control mode memory,” i.e. for keeping track of the control 
modes entered prior to the current time in the system ’s evolution. 
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Each set V (e INIT ' ,a , e DE5T;a ) can be computed u sing the Depth-First 


1990 


Search algorithm, whose computational cost is known flCormen et ah 
section 23.3) to be 0(|V| + |i?|)0. By using parallel computation in the 
Depth-First Search, however, one may reduce the physical running time by a 
factor proportional to the number of processors available. Another reduction 
of computational cost of routing in the context of the SSRP is attempted 
in f Sadovskv . to appear]') bv impo s ing a weighted graph structure on the 
Cartesian product flKorn and Korn . 196ll . 12.7-1) 


H V (e INIT ' a , 

aeA 


e DEST-,a 


in such a way that the shortest paths give the routings for the most desirable 
control modes p . 

These considerations suggest the following stepwise procedure for finding 
an optimal solution to an SSRP, in which the worst case-of having to enu- 
merate all the control modes explicitly-is assumed (step [2] in the procedure). 
Upper bounds on the computational cost are provided, in []’s, where possible. 

1. Compute all the sets V ( e INIT '^ a ^ e DEST '^_ [if this is accomplished by 
running Depth-First Search for each aircraft a G A, hence the compu- 
tational cost is 0(|*A|(|U| + |£'|)). The physical running time, however, 
can be reduced by using parallel processing; the reduction would be by 
a factor proportional to the number of processors.] 

2. Enumerate all the control modes p. [An upper bound exponential in 
the 1*4.1 .] 

3. For each p, compute an optimal solution to the corresponding Stacked 
Scheduled Routing Problem, and the cost of that solution. [The 


7 Here V and E denote, respectively, the vertex set and the edge set of the multigraph 
that models the airspace, as defined in item [J in the numbered list o f sectionU The symbol 
O has the following meaning (see, e.g.. ICormen et al~ ( 199Clh and ( Korn and Kornl [l96ll 
section 4.4-3)): if f(n) and g{n) are functions of the positive integers n, and if starting 
with some integer no one has 


|/(n)| < (some constant)g(n) for all n > no, 

then one calls g(n) an upper bound for f(n ) as n — >• +oo and writes “/(n) = 0(g(n)) as 
n — > +oo,” or simply u f(n) = 0(g(n)).” 
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computational cost depends on th e part icular choice of the computa- 
tional method; see, for example, Pytlak (119991 1 . The numerical re- 
sults in this paper were produced using a solver based on Sequential 
Quadratic Programming (SQP); see section I5l3l Since SQP is iterative, 
the task of estimating computational costs runs into the difficulty of 
estimating the rate of convergence, which depends on the specific initial 
guess and gradients arising in the execution of the method.] 

4. Select a control mode p* such that C^* < C ' M for all //, and declare 
the corresponding optimal solution the optimal solution for the SSRP 
(definition 14. 2p . [0(the number of control modes).] 


5. Sample numerical computations for the Stacked Scheduled 
Routing Problem 

5.1. Assumptions and notational conventions 

All the route networks appearing in the numerical examples of this sec- 
tion are graphs; namely, no two vertices are connected by more than one 
edge. This allows to specify each path as a sequence of vertices, rather 
than of edges. To simplify the computations, the separation requirements 
for each pair of aircraft are assumed symmetric. The numerical code admits 
a straightforward, albeit somewhat cumbersome algebraically, generalization 
that will dispense with this assumption. 

In panels (a, b) of Figures |9j [hQ fl2lfl7l the computed trajectories and 
controls are plotted using the symbols described in Table [2] In panels (A) 
(and, when present, (D)) of these Figures, as well as in Figures 0 [TQ1 the 
route networks are depicted as directed graphs; aircraft, as points labeled 
with values of a, each point serving as a center of a circle with radius equal 
to half the required pairwise separation. In all plots, axis label p refers, in 
agreement with the ab ove, to physical time. All plots were generated using 
the Matlab software MATLAbI ( 2010 ). 

The first two of the examples in section HTT1 are “abstract,” in the sense 
that no particular application is specified for them. Thus, the units of length 
and physical time are left unspecified. Application and units are, however, 
specified for the third example. The required destination y DEST ’ a is in each 
case the end of the path assigned to aircraft a in control mode p: 


V 


DEST;ct _ 


= 0 . 
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Quantities Plot symbol 


y\:*\ 

y 2 q ,* 2 q + 

y\,s\ 


> for q — 1 , . . . , A. 


Table 2: The legend used in Figures l9l ITT1 fl2lfT6l below. 


5.2. The cost function 

The cost function used in the following examples is ([TUP which, on convert- 
ing the scheduled routing problem to a Stacked Scheduled Routing Problem 
f definition 14. 2p . takes the form 


E 




a , 

Sq 


Pqdr. 


<7=0 q '> q + 1 ' 


5.3. Computational methods for finding optimal control strategies 

In each control mode, the corresponding optimal control problem can be 
approached using a numerical method from any suitable family; see, e.g., 
Betts ( 1998 1 for a survey of such methods and Pvtlakl ( 1999h for a thorough 
exposition. For reasons of convenience at the time of this research (easy ac- 
cess to code which is open-source, platform- independent, and self-sufficient), 
the authors used the OCP solve r (IFabie n. 20 1 If) . which is based on Sequential 
Quadratic Programming ( SQP ) Pvt lakl (119991) . 


5.f. Numerical examples 

5 . 4 .I. Three aircraft, one control mode 

The route network for this example is the 2-dimensional directed graph 
shown, together with the initial locations of the three aircraft, in Figure |8l 
Each aircraft can traverse only the edge on which it is positioned initially, 
hence only one control mode p arises. 

The speed ranges are given by 

a 12 3 

s MlN-,a 0.3 0.3 0.4 

S MAA> L5 Q_g 0 # 9 
A 6 
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Initial state 



Figure 8: The route network and the initial locations of the aircraft in the example of 
section 15.4.11 


The required times of arrival at destination are 

ct 12 3 

t VEST ' a 2.0 3.0 4.0 

The minimal required separation is 0.3. 

The numerical solution computed for the control mode /i is shown in 
Figure El 

Remark 5.1. The numerical data for this example show that an optimal 
control strategy (speed profile) is piecewise constant. For solutions in which 
no two aircraft are ever at minimal separation , this structure of an optimal 
control strategy is consistent with the Pontryagin Maximum Principle. 

5.4-2. Two aircraft, four control modes, wide speed ranges 

The route network for this example is the 2-dimensional directed graph 
shown, together with the initial locations of the two aircraft, in Figure [TOj 
The two paths considered here are 


pi : v 1 ,V2,v 8 ,v 9 ,v 10 ,v 6 ; p 2 : vi,v 2 ,v 3 ,V4,v 5 ,v 6 


The speed ranges are given by 


a 

MlN-,a 

MAX-, a 
b V 


1 2 
0.6 0.6 
1.4 1.4 
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Figure 9: Numerical solutions (for the only possible control mode) in the example of 
section 15.4.11 (A) State trajectory vs. time. (B) Control strategy vs. time. (C) The 
positions of the aircraft in the route network at p q {l) for q = 1. (D) The positions of 
the aircraft in the route network at p q ( 1) for q = 2 (aircraft a = 2 has just arrived at its 
destination and exited the system). 


The required times of arrival at destination are 

a 1 2 

t DKST,a 28.3 36.1 
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Initial state 



-fa -5 0 5 10 

Figure 10: The route network and the initial locations of the two aircraft in the example 
of section 15.4.21 

The minimal required separation is 1.0. 

The control modes are as follows. 


a 

Pi (a) 

p 2 (a) 

p 3 (a) 

p 4 (a) 

1 

Pi 

P2 

P2 

Pi 

2 

P2 

Pi 

P2 

Pi 


The respective numerical solutions for these modes are shown in Figures 

EH 

The total costs for each control mode are as follows: 


control mode 

Pi 

p2 

P3 

P4 

total cost 

37.5823 

38.8925 

45.4652 

31.5979 


Therefore, the optimal routing and control strategy are achieved in control 
mode /i 4 . 

Remark 15.11 applies to the numerical results shown for this example. 

5. 4-3. Two aircraft, four control modes, narrow speed ranges 

In this example, the route network, 3-dimensional, models a terminal 
airspace, which consists of two arrival paths merging into the same final 
approach segment. The route network and the initial positions of the two 
aircraft is shown in Figure El Here the edge (u 6 . vf) is the final approach to 
a runway. For realism, one unit of arc length is taken here to be 3 nautical 
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(A) (B) 


End of time interval #1 : p = 28.3 



(C) 

Figure 11: Numerical solutions for control mode pi in the example of section 15.4.21 (A) 
State trajectory vs. time. (B) Control strategy vs. time. (C) The positions of the aircraft 
in the route network at p q { 1) for q = 1. 


miles, and the unit of speed is taken to be 200 knots, a typical speed allowed 
in U.S. terminal airspaces at altitudes below 10,000 feet. 

The two paths considered here are 


pi : V \ , v 2 , v 8 , v 9 , vio, v 6 ; p 2 : Vi, v 2 , v 3 , n 4 , v 5 , v 6 


28 




(A) (B) 

End of time interval #1 : p = 28.3 



-fa -5 0 5 10 

(C) 


Figure 12: Numerical solutions for control mode p 2 in the example of section 15.4.21 (A) 
State trajectory vs. time. (B) Control strategy vs. time. (C) The positions of the aircraft 
in the route network at p q { 1) for q = 1. 


The control modes are as follows. 


a 

Pi(a) 

P 2 («) 

H 3 (a) 

Pi (a) 

1 

Pi 

P2 

P2 

Pi 

2 

P2 

Pi 

P2 

Pi 
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(A) (B) 



(C) 

Figure 13: Numerical solutions for control mode p 3 in the example of section 15.4.21 (A) 
State trajectory vs. time. (B) Control strategy vs. time. (C) The positions of the aircraft 
in the route network at p q { 1 ) for q = 1 . 


The speed ranges are given by 

ct 1 2 

0.8 (= 160 kts) 0.8 
s“- Y; " 1.2 (= 240 kts) 1.2 
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(A) (B) 


End of time interval #1 : p = 28.3 



Figure 14: Numerical solutions for control mode p 4 in the example of section 15.4.21 (A) 
State trajectory vs. time. (B) Control strategy vs. time. (C) The positions of the aircraft 
in the route network at p q { 1) for q = 1. 

The required times of arrival at destination are 

ct 1 2 

t uEST-, a 21.2 (= 0.03 hrs) 28.1 (= 0.04 hrs) 

The minimal required separation is 1.0 (= 3 NM). 
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Initial state 



Figure 15: The 3-dimensional route network and the initial positions of the aircraft 

(centers of the circles in horizontal planes) , for the numerical example of section 15.4.31 
Here the route network is a terminal airspace. The radius of each circle is half the minimal 
required pairwise separation of 3 NM. 


In control modes /i 2 ,/i 3 , the control problem has no solution, since the 
path lengths are such that the imposed speed ranges prevent at least one 
aircraft from reaching its destination at time. The numerical solution for 
control mode /ii is shown in Figure fl6l and incurs a total cost of 49.2396. In 
/i 4 , the computed optimal control strategies prescribe the constant speeds 

^ = 1 . 03 , < = 0.82 

and the solution (Figure fTTll incurs a total cost of 41.5594. Therefore, the 
optimal routing and control strategy are achieved in control mode ^ 4 . 
Remark 15. II applies to the numerical results shown for this example. 


6. Discussion 

The above modeling framework addresses the problem of navigating a set 
of aircraft in a route network, with constraints on the initial locations of the 
aircraft, required destinations, and required times of arrival at destination, as 
well as on the minimal pairwise distances between aircraft. We now discuss 
several directions in which the above model can be varied and generalized. 


6.1. A model that includes inertia 

Inertia can be included by treating both the y“’s and the s“’ s as state vari- 
ables, and the accelerations a“ as the control variables. The corresponding 
new state equations would assume the form 


Vu 




a. £ 
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(A) (B) 

End of time interval #1 : p = 21 .2 



(C) 


Figure 16: Numerical solutions for control mode in the example of section 15.4.31 (A) 
State trajectory vs. time. (B) Control strategy vs. time. (C) The positions of the aircraft 
in the route network at p q (l) for Q = 1- 


The rescaling of physical time to normalized time r G [0,1] formulated in sec- 
tion [472] would be carried out analogously. The resulting problem falls in the 
class of the kinodynami c m oti on planning problems; some related theoretical 
results can be found in .Tung and Christ (2008) and in references therein. 
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(A) (B) 

End of time interval #1 : p = 21 .2 



(C) 


Figure 17: Numerical solutions for control mode p.4 in the example of section 15.4.31 (A) 
State trajectory vs. time. (B) Control strategy vs. time. (C) The positions of the aircraft 
in the route network at p q { 1) for q = 1. 


6.2. Different choices of the cost function 

Of the vast number of possible choices of cost function, we briefly discuss 
two. 


• In situations where inertia cannot be neglected and acceleration must 
be the control variable (or among the control variables), it may be de- 
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sirable to keep the movement as smooth as possible, e.g. for passenger 
comfort or cargo safety. A cost function that would serve this goal is 
the sum of the integrals of the squared accelerations: 



• In situations where it is undesirable to impose times t DEST ’ a of arrival at 
destination as rigid constraints, and preferable to minimize the absolute 
differences between the actual arrival times f AAT ’ a , This goal would be 
served by taking the sum 



DEST,a 




J.AAT ;q 


) 


2 


as the cost. 


6.3. Asymmetric and anisotropic pairwise separation requirements 

In some applications, the pairwise separation requirements for aircraft can 
be asymmetric (i.e., dependent on the relative position of the two aircraft) or 
anisotropic (i.e., dependent on the specific location in the airspace), or both. 
An example of asymmetry is as follows. If two aircraft are consecutively 
in-trail and a t the same altitude, and th e leader’s and follower’s respective 
weight classes [Order JO 7110. 65U (2012) are Heavy and Small, the required 
separation is considerably larger than if the two aircraft types were in the 
opposite order. An example of anisotropy is the requirement of vertical sep- 
aration between two aircraft, where aircraft are required to maintain either 
1000 ft vertical separation or the prescribed lateral separation previously 
discussed; the resultant shape of an aircraft’s safety envelope is a cylinder. 
A mathematical form for such an anisotropic constraint would use, not the 
Euclidean norm, but one of the following form: putting a fc = (a£, a|, af.), 
the norm in the left-hand side of ([5D would be replaced by a “mixed” norm, 
Euclidean in the xy-plane, and the max-norm along the height z. If r were 
the minimal horizontal separation, imposed when the two aircrafts’ altitudes 
differ by less than a quantity h , the mixed norm and the corresponding sep- 
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aration constraint would be 


Cia i c ' 2 a 2 1 1 mixed 


:= max (a af - c 2 af) 2 + (ad[ - c 2 a y 2 f , \ |ciaf - c 2 a|| j 

> 1 

The computations in section [3] would have to be modified accordingly and 
would no longer have the closed quadratic form. 
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